The Landau-Lifshitz/Looyenga dielectric mixture expression and its self-similar fractal nature 
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In this paper, dielectric permittivity of dielectric mixtures is discussed in view of the spectral density rep- 
resentation method. A distinct representation is derived for predicting the dielectric properties, permittivi- 
ties e, of mixtures. The peculiar presentation is based on the scaled permittivity £ = (e e — e m )(ei — £m) , 
where the subscripts 'e', 'm' and 'i' denote the dielectric permittivities of the effective, matrix and inclu- 
sion media, respectively [Tuncer E 2005 J. Phys. Condens. Matter 17 LI 25]. This novel form of represen- 
tation is the same as the distribution of relaxation times formalism in dielectric relaxation. Consequently, 
we propose an expression for the scaled permittivity, which is the same as one of the extensively used di- 
electric dispersion expressions, known as the Havriliak-Negami empirical formula. The scaled permittivity 
representation has potential to be improved and to be implemented in to the existing analyzing routines for 
dielectric relaxation data to extract the topological/morphological description in mixtures. In order to il- 
lustrate the strength of the representation and confirm the proposed hypothesis, Landau-Lifshitz/Looyenga 
expression is selected, and the structural information of the mixture is extracted. Both a recently devel- 
oped numerical method to solve inverse integral transforms and the proposed empirical scaled permittivity 
expression are employed to estimate the spectral density function of the Landau-Lifshitz/Looyenga ex- 
pression. In the simulations the concentration q of the inclusions phase are varied. The estimated spectral 
functions for the mixtures with different inclusion concentration compositions show similar spectral den- 
sity functions, composed of couple of bell-shaped distributions, with coinciding peak locations. We think 
therefore that the coincidence is an absolute illustration of a self-similar fractal nature of mixture topology 
(structure) for the considered Landau-Lifshitz/Looyenga expression. Consequently, the spectra are not al- 
tered significantly with increased filler concentration level-exhibit a self-similar spectral density functions 
for different concentration levels. Last but not least, the calculated percolation strengths also confirm the 
fractal nature of the systems characterized by the Landau-Lifshitz/Looyenga mixture expression. We con- 
clude that the Landau-Lifshitz/Looyenga expression is therefore suitable for complex composite systems 
that have hierarchical order in their structure, which confirms the finding in the literature. 

PACS numbers: 77.22.-d, 78.20.-e, 77.22.Ch, 77.84.Lf, 02.70.Hm, 02.70.Uu, 05.45.Df, 07.05.Kf, 61.43.-j 
Keywords: Dielectric properties of solids and liquids , Optical properties of bulk materials and thin films , Permittivity 
(dielectric function) , Composite materials , Spectral methods , Applications of Monte Carlo methods , Fractals , Data analysis: 
algorithms and implementation; data management , Disordered solids 



INTRODUCTION 



Electrical properties of composite materials have at- 
tracted researchers to seek a relation between overall com- 
posite properties and intrinsic properties of the parts form- 
ing the mixture (constituents) and their spatial arrangement 
inside the mixture Q S S B S IE S IE M EO. Ell Q E3l 
Mixture formulas based on analytical and effective medium 
approaches were developed, such that for various arrange- 
ment of inclusions predi cting th e dielectric properties of 
composites was plausiblel 14i ll5L fT^. fT7l fT^l . A deep un- 
derstanding of dielectric mixtures would be of great value 
(i) to be able to calculate either the dielectric constant of a 
mixture of substances of known dielectric constants or, (ii) 
knowing the dielectric constants of a mixture of two com- 
ponents and that of one of the components, to calculate 
the dielectric constant of the other 1 12], or (iii) even know- 
ing the dielectric constants of a mixture of two components 
and that of two components to estimate the morphology of 
the mixture ll^l2fjLl2lll . In late 1970's, Bergman cleverly 
showed that one can separate the geometrical contributions 
from the pure dielectric response of a composite if and only 
if the dielectric properties of the constituents were known 



Milton corrected errors in the Bergman's orig- 
inal derivation 1I25H26II27I1 and later Golden and Papanico- 
laou I2ll29ll gave the rigorous derivation for the spectral 
representation theory. Recently, the present author has il- 
lustrated similarities between the dielectric relaxation and 
dielectric response of dielectric mixtures using the spectral 
density representation, the origin of similarities is very sig- 
nificant to comprehend physics of dielectrics lEoll . 

The concept of having the cognition of the structure of 
composites, how the phases are arranged, is very useful in 
materials design, because special materials can be manu- 
factured with the knowledge of structure-property relation- 
ship. For regular arrangement of phases, there exists equa- 
tions based on theoretical calculations on simple enough 
geometries. However, fractal structures are abundant in na- 
ture 130I1 . therefore to comprehend the materials properties 
with fractal structure has been a challenge for researchers 
for some decades. The fractal geometry or systems indi- 
cating hierarchical order has been one of the interesting 
topics in applied and theoretical (mathematical) physics 
l3fjl I3II I32I l33i 13]. As an example, electrical proper- 
ties of metal aggre gate s in insultin g m atrix media were 
studied extensively M HE El MM EE EH EI- 
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In these studies either the dimension of the electrical net- 
work or the system, or the resonance frequency of the elec- 
trical impedance was used as a measure to indicate the 
fractal dimensions. When real systems are taken into ac- 
count, the structural information is on the other hand usu- 
ally obtained by optical/microscopic techniques, which are 
later analyzed to estimate the fractal dimensions. One can 
as well utilize a mixture formula to model the electrical 
properties of the composite system in hand, such that the 
model contain structural information, e.g. there exist ef- 
fective medium theories for composites with spherical and 
ellipsoidal inclusions J2, 13, 145II . In the present paper, 
we employ the spectral density representation, which is a 
general representation for composites, to resolve the geo- 
metrical description of a model system described by the 
Landau-Lifshitz/Looyenga (LLL) effective medium for- 
mula [461 14711 . or in other words we challenge the physical 
significance of the LLL expression. The LLL expression 
was used to describe the dielectric properties of dispersive 
systems composed of powders or exhibiting porous struc- 
ture llliilSE3EilE3JllElElESEiElElE3 

It was even showed |4<|] that the LLL formula was 
more reliable when mixtures contained strongly dissipa- 
tive particles and compared to others like Maxwell Gar- 
nett (MG) |( ffij < 53j)- Bruggema n I64II. e tc. (see for example 
Refs. flllBlIHElllSElilinnil for other formulas). 

In this paper, we first present that the spectral density 
representation can in fact be written in a novel, more ele- 
gant, form that can be implemented in already existing di- 
electric data analysis techniques Later 
we use the presented novel notation and the numerical pro- 
cedure to solve inverse problems I19l l2lL l7(A \l ill to test 
our hypothesis. We proceed to achieve our goal by consid- 
ering the LLL ll46l l47ll expression for dielectric mixtures 
and discuss the significance of the presented approach on 
the mixture expression. 

The paper is organized as follows, first we present the 
spectral density representation for a binary mixture in % 
in this section the similarities between the dielectric re- 
laxation in dielectrics and dielectric permittivity of binary 
mixtures are illustrated. % describes the dielectric data 
representation and gives hints for analyzing impedance 
data of mixtures. The numerical method to solve the in- 
verse integral is also presented explicitly for the interested 
readers in % The numerical data generation and Landau- 
Lifshitz/Looyenga expression are presented in % The com- 
parison of the results obtained by the inverse integral so- 
lution and the proposed conventional dielectric dispersion 
expression are given in ^ Conclusions are in ^ 



SPECTRAL DENSITY REPRESENTATION 

In the spectral density representation analysis of binary 
mixtures, the dielectric permittivity of a heterogeneous (ef- 
fective) medium, is expressed as l22lE3il24lE5[E5E7lE8i 



H|7l|7l|7l, 



1 + . 



^ m 



1 + 



G(x)dx 



7^-1 



+ X 



where e e , e m and e-, are the complex dielectric permittivity 
of the effective, matrix and inclusion media, respectively; 
q and x are the concentration of inclusions and the spec- 
tral parameter, respectively. The function G(x) is the spec- 
tral density function (SDF), possesses information about 
the topological description of the mixture. Eq. can be 
arranged in a more elegant form, cf. § D as follows, 



G(x)dx 



(2) 



10 l + Airn^ 

where, A im = e-, — e m , and £ is the complex and frequency 
dependent 'scaled' permittivity, 

S & Sjnn 



(3) 



The constant £ s in Eq. |2] is complex and depends on the 
concentration and structure of the composite, its real part 
is related to the so called the 'percolation strength' l72i 
l75ll . The mathematical properties and conditions that SDF 
satisfies are presented in ^I24ll72ll73ll74ll76ll. 

Eq. (|3j is a very similar expression to the distribution of 
relaxation times (DRT) representation of a broad dielectric 

dispersion (relaxation) Efl E3, El El El HH EH US HI 

f°° G(r)dr 
e{uj) = £oo + Ae - \ ' (4) 

where e, and Ae are the complex permittivity, permit- 
tivity at optical frequencies and dielectric strength, respec- 
tively; and 1 = The quantities cj and r are the angu- 
lar frequency and relaxation time, respectively. The distri- 
bution function for the relaxation times is G(r). Compar- 
ison of Eq. (|2j and (@J demonstrate that both the DRT and 
the scaled permittivity of SDF are actually the same. How- 
ever, a new complex parameter w in SDF w = e~ 1 A im 
corresponds to the pure complex frequency iuj in DRT rep- 
resentation, and the real number constant is a complex 
number £ s in SDF representation. In addition the spec- 
tral parameter x corresponds to the relaxation time r in the 
DRT. Finally, the dielectric strength Ae in the DRT repre- 
sentation is related to the concentration of inclusions in the 
SDF. 

Due to the presented similarities or in other words the 
analogy, methods developed for dielectric data analysis ll65l 
I68L l69ll can be applied to the scaled permittivity £ of SDF 
lll9l l20t I2II1 . For example, one of the most employed 
dielectric dispersion expressions, known as Havriliak- 
Negami empirical expression [85], can be used to analyze 
the scaled complex dielectric permittivity data of a mixture, 



i(w) = 6 + 



[1 + {wx) a \ 
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FIG. 1: The Argand diagram of resistivity for a MG compos- 
ite in log-log scale, the effective permittivity is calculated with 
e e = e^ G (w;e m , e h 0.3,3), withe™ = 2 + lO -12 ^^) -1 and 
e\ = 10 + 10 -10 (k w) -1 . The arrows indicate the direction of 
increasing angular frequency uj. The insets (a) and (b) are the 
real and the imaginary parts of permittivities; e e , e m and e, are 

the permittivities for the matrix ( ), inclusion ( ) and 

effective ( ) media. The inset (c) is the Argand diagram of 

resistivity in linear-linear scale. Logarithmic scale is base 10. 



where a and are parameters of a general distribution 
function and w is the scaled complex frequency. 

When a = = 1, a Debye-type relaxation is observed 
in the dielectric dispersion representation 1 86] . In the case 
of spectral density representation, however, the Maxwell 
Garnett approximation is obtained for a = = 1 and 
x = ( 1 — q) Id, since G(x) = 5[x — (1 — q)/d] 

IlllH EZal, where d is the dimensionality of the 

system. 

One can therefore in principle write a new more general 
empirical mixture formula [87] by isolating the dielectric 
permittivity e e of the composite in Eq. (|3j and substituting 
it in Eq. (|5J, 



e e = e m + A im <^ i s + 



[1 + {wx) a Y 



(6) 



Finally, note that the second, fractional expression inside 
the curly parenthesis in Eq. © can be exchanged with any 
one of the dielectric dispersion relations existing in the lit- 
erature SlEimmmiiiJIiillii. 



REPRESENTATION OF DIELECTRIC DATA 

The dielectric function for a <i-dimensional (or com- 
posite with arbitrary shaped inclusions) is defined as fol- 
lows with Maxwell Garnett (MG) expression for a com- 



-0.5 



o 



-1.5 













* 

/ 


\\ 






w 

t 

t 


\ \ 


\ 


t 

i 


t 




\ 

\ 

\ 


t 

1 

1 






\ 

\ (a? 



-2-1 1 2 

log(co [rads -1 ]) 



0.4 



I 0.2 













(b) 



























0.2 0.4 



0.6 
X'(0>) 



0.8 



FIG. 2: (a) Dielectric susceptibility \ as a function of angular fre- 
quency uj; the real and the imaginary parts are presented with the 

solid ( ) and dashed ( ) lines, respectively. The arrow 

indicates the direction of increasing angular frequency uj. (b) The 
Argand diagram (Cole-Cole plot) of susceptibility. Logarithmic 
scale is base 10. 



posite IHIIHI 

ef G (w; e m , e h q, d) = e m + 



e m dqA in 



[l-q)A- m + de„ 



(7) 



The dielectric data of the composite can be expressed in 
one of the four immittance representations [6g, |68j, 17811 . 
(z) the complex resistivity p(u>); (ii) the complex mod- 
ulus M(u) = iu)Eqp{uj); (Hi) the complex permittiv- 
ity e = [M(lj)] _1 ; and (iv) the complex conductivity 
ct(lu) = iue e{uj) = [p^)]^ 1 . When we are deal- 
ing with frequency dependent dielectric properties of com- 
posites the effective conductivity of the composite can 
sometimes influence the imaginary part of the dielectric 
function-hinders the dielectric losses due to the interfa- 
cial polarization-as shown in inset (b) of Fig. ^ In suc h 
cases it is more appropriate to use the complex resistiv- 
ity representation (plot) as shown in the Argand diagrams 
in log-log and linear scales in Fig. Q and respectively. 
The ohmic conductivity cr e ( j c (or resistivity /9 e ,d c ) of the 
material can be estimated from the complex resistivity Ar- 
gand plot. However once the conductivity contributions 
are cleared from the immittance or dielectric data, us- 
ing the estimated resistivity value in Fig. ^ as w ~~ * 0, 
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the pure dielectric dispersion (permittivity) would be ob- 
tained. In addition the high frequency dielectric permit- 
tivity e(u — > oo) = can be further subtracted from 
the data to obtain the pure dielectric polarization (suscep- 
tibility x) °f the composite as presented in Fig. |3 X = 
X' — ix" = £ — Soo + kt^Equ) -1 . The imaginary part of 
X has a peak around to ~ 1 + rads -1 . The linear scale plot 
of the susceptibility x is a sem i circular curve, cf. Fig.lJJ}, 
as seen in Fig.^. 

If we now consider the frequency dependent properties 
of the scaled permittivity £ for the considered MG com- 
posite, the real part and the imaginary parts are similar to 
the dielectric permittivity e e of the composite, cf. inset in 
Fig-El which shows the Argand diagram of the scaled per- 
mittivity £, observe that the increasing frequency is in the 
opposite direction when compared to the Argand plot of the 
susceptibility in Fig. 03. The real part of £ is a mirror image 
of 5R(e e ). Unlike the imaginary part of e e (due to the ohmic 
losses), shows a clear peak around u ~ 1~ rads -1 . 
The shift in the origin position in the Argand diagram in 
the inset of Fig.|3]is related to the percolation strength £ s , 
which is close to zero, and the concentration of the inclu- 
sions. There are similarities between the susceptibility and 
the scaled permittivity plots when the same frequency uj 
axis is used, however, note that the scaled frequency w for 
the scaled permittivity is a complex quantity. We therefore 
illustrate the dependence of the real angular frequency uj as 
a function of w in a 3D curve-plot in Fig.|4] In addition the 
real and the imaginary parts of the scaled permittivity £ are 
shown in Fig. |5] as a function of w. As shown in the fig- 
ure the actual dependence of £ on w is more complicated 
than e on uj. On the contrary, this dependence can be used 
to estimate the SDF for a given system, which is explicitly 
given with a numerical procedure in the next section. 



NUMERICAL ESTIMATION OF SPECTRAL DENSITY 
FUNCTION 

The derived spectral density expression in Eq. UJ is a 
Bolter equation [93], which is a special form of the Fred- 
holm integral equations iEmIi . Such equations are usu- 
ally considered to be ill-conditioned because of their non- 
unique solutions. However, the approach used here and re- 
cently presented elsewhere ll9U67ll70l 84] leads to unique 
solutions by means of a constrained least-squares fit and 
the Monte Carlo integration methods. Some other ap- 
proaches were also suggested to solve the spectral density 
function in the literature f^Fll7ll7ll7^l9lll9^l97[l9il. 
The presented numerical method to solve inverse inte- 
gral transforms has previously been used in different prob- 
lems El S O |70l Ell & In this particular ap- 
proach, the integral in Eq. (|3j is first written in a sum- 
mation form over some number of randomly selected and 
fixed x„-values, x„ 6 [0, 1], where n is less than the to- 
tal number M of experimental (known) data points in the 
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FIG. 3: The real 
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parts of scaled permittivity £ in the Maxwell Garnett approxima- 
tion. The inset is the Argand diagram for £, the arrow shows the 
direction for increasing angular frequency uj. The inclusions are 
spherical d = 3, and q — 0.3. Logarithmic scale is base 10. 




FIG. 4: Dependence of the real frequency u> as a function of com- 
plex 'scaled' frequency -cu in the spectral representation. Loga- 
rithmic scale is base 10. 



complex scaled permittivity, £, 



f=«.+£ T 



+ wx r , 



n < M 



(8) 



This converts the non-linear problem in hand to a linear 
one with g„ being the unknowns, weights of the randomly 
selected x n values. In the present notation g is q G. Later, 
a constrained least-squares algorithm is applied to get the 
corresponding g-values and £ s , 



mm 



^[e-£fl] 2 and g>0 



(9) 
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where & is the kernel-matrix, 



/ 1 K n K 12 
1 K 21 K 22 
1 K 31 K 32 

V: : : 



(10) 



Here Ky = [1 + WiXj]^ 1 , index z runs on the angular fre- 
quency points i = 1, . . . , M, and index j runs on the ran- 
domly selected x values, j = 1, . . . , n. The parameters £ 
and g in Eq. (|9j are column vectors, respectively, the scaled 
permittivity calculated from the experimental (known) data 
and the searched spectral density, 



/fMwi)]\ 



V 



and 



/ 



gi 

g2 
g3 

V : J 



(11) 



In our numerical procedure, we perform many minimiza- 
tion steps with fresh, new, sets of randomly selected Xj- 
values. The g., -values and £ s obtained are recorded in 
each step, which later build-up the spectral density dis- 
tribution g and a distribution for the percolation strength, 
£ s . For a large number of minimization loops, actu- 
ally the x-axis becomes continuous — the Monte Carlo 
inte grati on hypothesis — contrary to regularization meth- 
ods |82l EBl]. In the presented analysis below, the to- 
tal number of randomly selected x values are ~ 10 6 . The 
number M of data points are chosen 24, and the number 
N of unknown g-values are 22. 

Application of the numerical procedure to the Maxwell 
Garnett expression, impedance data of porous rock-brine 
mixture and two-dimensional 'ideal' structures has previ- 
ously been presented elsewhere [ 19, The estimated 
spectral density functions for the MG expression were delta 
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FIG. 6: The Argand diagram of resistivity for a LLL compos- 
ite in log-log scale, the effective permittivity is calculated by 
e e = £e LL (^;£m, Ej, 0.3, 3) with e m = 2 + 10" 12 (i£ ui) ~ 1 and 
£i = 10 + 10~ 10 (teocoy 1 . The arrows indicate the direction of 
increasing angular frequency w. The insets (a) and (b) are the 
real and the imaginary parts of permittivities; e e , £ m and e, are 

the permittivities for the matrix ( ), inclusion ( ) and 

effective ( ) media. The inset (c) is the Argand diagram of 



resistivity in linear-linear scale. Logarithmic scale is base 10. 



sequences 1 99] as expected, without any significant perco- 
lation component, because of the estimated concentration 
q, cf. Eq. (I20t . In the next section, we apply the numer- 
ical procedure to the Landau-Lifshitz/Looyenga Efl l47ll 
expression in order to better understand the nature of di- 
electric mixtures, which obey this relation. 
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FIG. 7: (a) Dielectric susceptibility x as a function of angular fre- 
quency lo; the real and the imaginary parts are presented with the 

solid ( ) and dashed ( ) lines, respectively. The arrow 

indicates the direction of increasing angular frequency u>. (b) The 
Argand diagram (Cole-Cole plot) of susceptibility. Logarithmic 
scale is base 10. 



LANDAU-LIFSHITZ/LOOYENGA EXPRESSION 

Landau and Lifshitz Er3l and Looyenga indepen- 
dently, using different approaches developed an expres- 
sion for dielectric mixtures, implies the additivity of cube 
roots of the permittivities of the mixture constituents when 
taken in p ro p ortion to their volume fractions (see Refs. 

EliiliiEflElESElElElEEElEiEllSEl for 

examples). 



[e e LL >;e m , e h q)] 1 ^ = (1 - q)e^ + qer (12) 

This expression is extensively used in the literature for 
powdered materials and optical properties of material mix- 
tures l ! 2i 14911 . In the following calculations, we choose the 
same values for the dielectric functions of the phases as be- 
fore, cf. Fig. n The concentration q of inclusion phase is 
varied between 0.1 and 0.9 in the simulations. 

The extracted spectral functions are like distributions, 
and they are analyzed by means of comparing th em with a 
known di stribution. We apply the Levy statistics lIlQOtllOlL 
ll 02t ll 03h . which is widely used for interacting systems in 
different research fields EH EH ES GUI EH 
Il07lll08| l. The Levy stable distribution is a natural general- 
ization (approximation) of the normal (Gaussian), Cauchy 
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FIG. 8: (a) The real and (b) the imaginary parts of the di- 
electric permittivity calculated with Eq. H2\ for q = 0.3, the 
permittivities of the matrix and the inclusion phases 
2 + 10~ 12 (i£ ^r 1 andei = 10+10 _lo (ie o w) _1 . (c) Argand di- 
agram of scaled permittivity without out the percolation strength 
contribution. The inset in (a) shows the dielectric susceptibility 
after the subtraction of ohmic conductivity and the permittivity 
at high frequencies. The inset in (c) shows the Argand plot of £ s 
for each point, the values are very narrowly distributed. The lines 
with points are the simulated data of Eq. i\2l and the symbols (o) 
are the dielectric response calculated with the estimated spectral 
density function from the proposed numerical algorithm. There 
is a very good agreement between simulated and analyzed data 
sets. Logarithmic scale is base 10. 



or Lorenz and Gamma distributions. It is used when ana- 
lyzing sums of independent identically distributed random 
variables by a diverging variance. Its characteristic func- 
tion is expressed as 

L(x;2l,/i, 7 ,C) =2l|exp{-|C(x- / u)r}| (13) 

Here, 7 is the characteristic exponent (7 > 0), /i is 
the localization parameter, £ is the scale parameter and 
21 is the amplitude. The special forms of Eq. d!3t are 
the Gaussian [L(x; A, /x, 2, £)], the Lorentz or Cauchy 
[L(x; A, [i, 1, ()] and Gamma [L(x; A, /j,, 1/2, ()] dis- 
tributions. Different forms of probability density functions 
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FIG. 9: Spectral density functions of Landau-Lifshitz/Looyenga 
equation for two concentrations; (a) q — 0.2 and (b) q — 0.3. 
The solid lines ( ) are Levy distributions adapted to the esti- 
mated peaks. The numerical data is presented with symbol, error 
bars. The numbers on the peaks denote the significant peaks. 
Logarithmic scale is base 10. 



for Levy statistics exists, we have adopted a stable distri- 
bution used in the literature llOOLIlOlLlloll . We omit the 
imaginary parts in the characteristic function because of 
their insignificance in the results. 

In Fig. [5] the simulated dielectric permittivity with 
Eq. (I12t is presented in the complex resistivity level. The 
actual dielectric data is shown in the insets Fig. [(J and|5p. 
Compared to the MG expression in Eq. @, shown in Fig.^ 
the LLL expression does not show a knee-point as the com- 
plex resistivity decreases with increasing frequency. In ad- 
dition the linear scale Argand diagram of the LLL expres- 
sion in Fig.|6j; is not a perfect semi circle as the MG one. 
The dielectric susceptibility after the subtraction of the fre- 
quency independent parameters and a are presented in 
Fig .0 The losses, x", are non- symmetrical for the LLL ex- 
pression, which indicates that the actual dielectric response 
can be modeled by a non-Debye dielectric dispersion, e.g. 
the Havriliak-Negami expression, this is also visible in the 
Argand diagram of the susceptibility, Fig. 03. 




log (x) 



FIG. 10: Spectral density functions of Landau-Lifshitz/Looyenga 
equation for nine different concentrations, q — {0.1, . . . , 0.9}. 
The data are shifted for better comparison. There are six peaks 
resolved between x — 10~ 3 and x = 1. Interestingly, the posi- 
tions of the most probable spectral parameters of the six are not 
altered with increasing concentration. Logarithmic scale is base 
10. 



The reconstructed complex dielectric permittivity and 
the scaled permittivity are presented in Fig.|8j It is impor- 
tant to mention that the fitting is performed in the scaled 
permittivity level. While numerically calculating the spec- 
tral function the randomly selected spectral parameters x 
are picked between 10~ 3 and 1 in logarithmic scale; the 
pre-distribution is log-linear (for details see I70B ). In such 
an integration limit in Eq. I0, the distribution or the func- 
tional contributions of x-values lower than 10~ 3 are in- 
cluded in the percolation strength £ s . 

The estimated spectral density functions for two concen- 
tration levels, q = 0.2 and q = 0.3, are presented in Fig. [9] 
There are six visible peaks, which are labeled with numbers 
one to six from left to right, respectively. We do not take 
into account peak (zero) in the analysis. The estimated in- 
tegral of the distributions are presented on the graphs. Each 
peak is analyzed by the Levy distribution as mentioned be- 
fore, and the solid lines ( ) show the individual dis- 
tributions, cf. Fig. [9] It is remarkable that the estimated 
distributions are similar in form but shifted a little bit up 
with increase in concentration q. This behavior is an illus- 
tration of self-similar fractal nature of the considered com- 
posite system in the LLL expression, in which the topo- 
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logical arrangement is not changing significantly with in- 
creased concentration. To support this argument/statement, 
the spectral density functions of mixtures described with 
the LLL expression for nine different concentration are 
shown in Fig. E3 The shifting is performed with con- 
stant steps in concentration for clarity, however as observed 
the spectral function amplitudes are not increasing propor- 
tional with increasing inclusion concentration, cf. Table U 
and cf. Fig. The difference between peak positions 
is constant in logarithmic scale \Pj — (Pj-i\ ~ 0.47, cf. 

Fig.ini 

In Table |2 not only the positions but also the shape pa- 
rameters of the Levy expression are listed for the six peaks 
resolved for each concentration. The peak positions [i are 
not varied with increasing concentration of the inclusions 
phase. The amplitude 21 changes with increasing concen- 
tration for each peak. The scale parameter £ and the char- 
acteristic exponent 7 indicate some relation to concentra- 
tion, the their exact relations are not sought. In Fig. \^\ 
we show the change in the amplitude for four of the peaks 
with increase in the concentration of the inclusions. It is re- 
markable that around q = 0.6 the behavior of the spectral 
density functions change, the increase in the amplitude of 
the peaks with increasing concentration starts to decrease 
for increasing concentration as q > 0.6. In the inset the 
ratio of the three selected peaks to peak 1 are shown in the 
figure, a similar activity is also observed in the ratio. Since 
the ratio between the amplitudes of the selected peaks does 
not indicate a simple linear relation to concentration q, the 
topological description of the system can not be qualita- 
tively investigated. However, as mentioned previously the 
location and form of the spectral density functions indicate 
a huge resemblance to each other that they are in fact re- 
lated to the self-similar hierarchical nature of the compos- 
ites expressed by LLL expression. 

Finally the statistical analysis of the concentration cal- 
culated from the integral of spectral density function [inre- 
gral expression on the right-hand side of Eq. 0], and the 
value of the percolation strength £ s at each Monte Carlo 
step are presented as summarized in Fig. El an d El The 
number distribution of integral of spectral density function 
calculated at each Monte Carlo cycle is not centered at the 
actual concentration q taken but deviated. The deviation 
indicates that there is a percolation path, network structure, 
cf. Fig. El As expected from the definition of the spectral 
density function, Eq. (I20t . addition of the estimated con- 
centration q and percolation strength £ s , presented with er- 
ror bars and open symbols (□) and (0)> respectively, yield 
very close numerical values as the actual concentration q. 
The value q is calculated from the distribution. It is striking 
that the applied numerical method is capable of estimating 
the concentration of the filler material when it is not known 
in advance. 
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FIG. 11: Spectral density function values at constant x, 
[log(x) = -2.5 (★); log(a;) = -2 (.); log(x) = -1.275 
(▼) and log (a:) = —0.675 (■)]. The inset illustrates the ra- 
tios for log^log^) = -2.5]}/log{g[log(z) = -2]} (□), 
log{g[log(a;) = -2.5]}/log{g[log( a; ) = -1.275]} (V) and 
log{g[log(s) = -2.5]}/Iog{gpog(x) = -0.675]} (o). Logarith- 
mic scale is base 10. 




FIG. 12: Peak positions, fi — log(:r). The solid line is a linear fit 
with n = 0.470Pj — 2.99, where j is the peak number in Tabled 



APPLICATION OF THE HAVRILIAK-NEGAMI 
EXPRESSION 

Previously, it has been stated that the scaled permittivity 
£ in Eq. (|3j can be expressed as in the convensional form as 
in the case of dielectric relaxation [20], Eq. (|5j. In order to 
show and verify this statement we apply a complex nonlin- 
ear least-squares curve fit algorithm to the scaled permit- 
tivity of the LLL expression, denoted as £ LLL below. The 
error in the curve fitting procedure is used to quantify the 
fitness of the model function in Eq. ©. The error is calcu- 



9 



TABLE I: Fit parameters of Levy distribution, Eq.^] for the seven peaks estimated by the numerical algorithm. 



Peak 


q 


21 


A* 


7 


c 


Peak 


q 


a 




7 


c 




0.1 


0.01 


-2.47 


1.58 


6.18 


2 


0.1 


0.06 


-2.01 


0.89 


15.84 




0.2 
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-2.49 
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6.57 
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-2.50 


1.94 


7.53 




0.5 


0.26 


-2.03 


1.30 


11.43 




0.6 


0.04 


-2.49 


2.06 


7.81 
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0.28 


-2.03 


1.33 


11.46 




0.7 


0.04 


-2.49 
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11.84 
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0.24 
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1.26 
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7.20 
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1.15 
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0.01 


-1.62 


2.97 


8.23 
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0.1 


0.23 


-1.22 


0.58 


35.98 




0.2 


0.02 


-1.62 


2.70 


8.26 




0.2 


0.54 


-1.24 


0.52 


49.50 




0.3 


0.03 


-1.61 


2.21 


9.34 




0.3 


0.82 


-1.24 


0.52 


53.04 




0.4 


0.03 


-1.61 


2.27 


9.10 




4 


1 04 


-1.24 


56 


49 16 
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0.04 


-1.61 


2.51 


8.84 




0.5 


1.20 


-1.24 


0.57 


48.10 




0.6 


0.04 


-1.60 


2.22 


9.05 




0.6 


1.30 


-1.24 


0.60 


47.57 




0.7 


0.03 


-1.61 


2.58 


8.91 




0.7 


1.26 


-1.24 


0.60 


48.15 




0.8 


0.02 


-1.61 


2.34 


8.69 




0.8 


1.05 


-1.24 


0.62 


45.98 
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0.01 


-1.62 


2.43 


8.32 
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0.58 
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0.62 
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0.38 


-0.14 


0.81 


10.58 




0.4 


0.62 


-0.62 


1.72 


10.89 




0.4 


0.40 


-0.14 


0.80 


9.84 




0.5 


0.69 


-0.61 


1.58 


10.90 




0.5 


0.37 


-0.14 


0.97 


9.08 
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0.65 


-0.61 


1.63 


10.25 




0.6 


0.33 


-0.14 


1.02 


8.87 
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0.56 


-0.61 


1.66 


9.57 
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0.28 


-0.14 


1.08 


8.93 




0.8 


0.39 


-0.60 


1.80 


8.46 




0.8 


0.18 


-0.14 


1.33 


8.22 




0.9 


0.26 


-0.60 


1.50 


9.33 




0.9 


0.10 


-0.14 


1.34 


8.45 



lated as the sum of relative error at each point as follows, 



sft(£ LLL ) - 5ft(£ ) 



LLL"i 



+ 



3(£ 



3(£ 



TABLE II: Fit parameters of Havriliak-Negami| 85] expression 
2 to scaled permittivity, Eq. The error <£ is calculated with 



Here £ HN is the model expression of Eq. Q. The fit re- 
sults are listed in Table UTI where the model values for the 
spectral parameter x, concentration q and the percolation 
strength £ s are presented with over-lines as in the previous 
section. 

The fit results are shown in Fig.^Jfor q = 0.3 case in 
the similar form as Fig. |8] There is actually no particular 
differences between the two methods, except that the nu- 
merical techniques based on the Monte Carlo algorithm is 
capable of resolving individual peaks, as shown in the com- 
parison graph in Fig. ^] It is not clear for example from 
the Havriliak-Negami approach that the system indicate a 
self-similar fractal-like structure. One should note that the 
Havriliak-Negami distribution, cf. § B |67], presented below, 
even spread over spectral parameter values larger than one, 
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X 


q 


Ts 




0.1 


1.034 


0.438 
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0.095 


0.000 


1.19x10" 
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0.555 
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0.009 


7.28x10" 


-6 


0.3 
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0.635 


0.464 


0.297 


0.032 


2.45x10" 


-4 


0.4 


0.850 


0.454 


0.522 


0.345 


0.066 


4.27x10" 


6 


0.5 


0.800 


0.470 


0.494 


0.391 


0.128 


2.06x10" 


-5 


0.6 


0.828 


0.422 


0.493 


0.397 


0.216 


4.08x10" 


-6 


0.7 


0.822 


0.409 


0.475 


0.371 


0.342 


5.20x10" 


-6 
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0.815 


0.399 


0.456 


0.300 


0.510 


6.41x10" 


-6 


0.9 


0.810 


0.390 


0.438 


0.179 


0.727 


7.60x10" 


-6 



x > 1, which is not possible in the spectral density repre- 
sentation. However, in order to analyze the data it is very 
convenient and trivial to implement in available curve fit- 
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FIG. 13: The number density distribution of estimated 
/ g(x)dx = q for two concentration q levels, q — {0.3, 0.4}. 
The dashed vertical lines are the expectation values of q, which 
are 0.255 ± 0.065 and 0.310 ± 0.096 for concentration levels 0.3 
and 0.4, respectively. Logarithmic scale is base 10. 




FIG. 14: The expectation values of constants, q and £ s , in the 
spectral density approach. The integral of g(x) yields the concen- 
tration of the inclusions q, denoted by (□) symbols. The expected 
value of percolation strength is £ denoted by (o) symbols. Both 
values are also presented with error bars. Addition of the two ex- 
pected values f s + q, denoted with (*) symbols, leads to the ac- 
tual concentration q as in Eq. <20t . which is denoted by the solid 
line ( ). In the inset, the number distributions of £ s , perco- 
lation strenght, estimated in each Monte Carlo cycle for q — 0.3 
and g = 0.4 are illustrated, respectively, the most expected £ s are 
0.041 ± 0.006 and 0.083 ± 0.015 for the two concentrations. 



ting programs. Similar to Fig.^J the estimated concentra- 
tion q and percolation strength £ s from the parametric anal- 
ysis satisfy the condition in Eq. i20t as shown in Fig.1171 
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FIG. 15: (a) The real and (b) the imaginary parts of the dielectric 
permittivity calculated with Eq. 1 12t for q = 0.3 and the modeled 
response obttained with the application of the Havriliak-Negami 
expression to the scaled permittivity. The permittivities of the 
matrix and the inclusion phases are the same as in Fig. [8] (c) 
Argand diagram of scaled permittivity without out the percola- 
tion strength contribution. The inset in (a) shows the dielectric 
susceptibility after the substruction of ohmic conductivity and 
the permittivity at high frequencies. The lines with points are 
the simulated data and the symbols (o) are the data estimates of 
the Havriliak-Negami expression, Eq. Q. There is a very good 
agreement between simulated and analyzed data sets. Logarith- 
mic scale is base 10. 



CONCLUSIONS 

In this paper, we first derived an expression for the 
dielectric mixtures that resembles the distribution of re- 
laxation times representation in dielectric relaxation phe- 
nomenon. In the derivation we used the spectral density 
representation. It is shown that the extisting knowledge on 
the dielectric relaxation theory can be applied to the di- 
electric properties of composites. In order to confirm the 
hypothesis, both a method similar to estimate the distri- 
bution of relaxation times and an extensively used empiri- 
cal formula to express dielectric relaxation are employed to 
estimate the spectral density functions of composites sim- 
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FIG. 16: Comparison of the distributions sfg obtained by the 

Havriliak-Negami ( ) and the numerical approach ( ). 

Observe that the Havriliak-Negami distribution also considers 
spectral parameter x values larger than 1, log (a:) > 0. Loga- 
rithmic scale is base 10. 
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FIG. 17: The most probable values of constants, q and £ s from 
the application of the Havriliak-Negami expression. The integral 
of g(x) yields the concentration of the inclusions q, denoted with 
(□) symbols. The percolation strength is f, denoted with (circ) 
sysmbols. The sum of the two constants, £ s + q, denoted with 
(*), leads to the actual concentration, which is denoted with the 
solid line ( ). 



ulated with the LLL expression. The numerical method 
based on the Monte Carlo technique estimated couple of 
peaks which did not change their location in the spectra 
with increased concentration of inclusions. This static be- 
havior of the spectra resembles that there exists a hierarchi- 
cal structural order in the composite, as a result, we infer 
that the LLL expression is proper for systems with self- 



similar fractal nature, such as composites with colloid ag- 
gregates and porous materials. We have explicitely shown 
why the LLL expression could be applied to describe the 
dielectric properties of powdered and porous systems. Last 
but not least, the findings are significant to confirm the 
structure of composite systems, whose dielectric permit- 
tivities are described with the LLL expression, in the liter- 
ature. 



DERIVATION OF THE SIMPLE FORM 

Eq. (0 is expanded as follows, 



qA 



+ 9 



/ £e — £n 
£ m 



G(x)dx 

v ' _! — a 5) 



Let £j — £j = Aij , then 



A— + ^ ^ G ( a; ) A imdx 



£-m Jo ~t~ ^im^ 

Now, multiply both sides with e m , 

^ qG(x)A im dx 



A em = q AA- m + 



o l + e m 1 A im x 



+ x 



(16) 



(17) 



Finally, let A em /A im = £ and qA = £ s , we obt ain E q. (|2j. 
The properties of G and A are such that lF7^.l73lll09ll . 



A+ f G{x)dx = 1, 
J{)+ 



a \ A 9(1-9) 



+ 



d 



(18) 



(19) 



Here, d is the dimension of the system. When we consider 
our new notation then, 



6 + / g(x)dx 



1 (1 
xg(x)dx = — 



(20) 
(21) 



0+ 



HAVRILIAK-NEGAMI DISTRIBUTION FUNCTION 

Havriliak and Negami |85] have combined the works of 
Cole and Cole [89]' and Davidson and Cole [88] and have 
expressed the dielectric dispersion with an asymmetric for- 
mula as presented in Eq. 0. When a = (3 = 1, Eq. Q 
becomes the simple Debye and Maxwell Garnett equations 
for dielectrics or dielectric mixtures, respectively. Other 
interesting cases are when a = 1; Davidson-Cole expres- 
sion and wh en = 1; Cole-Cole expression. Havriliak 
and Negami [110] have used the distribution of relaxation 
times as expressed by Davidson and Cole [88], and have 
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substituted sexp(±i7r) for is. The solution for the dis- 
tribution relaxation times and the spectral density function 
g(s), then, becomes, 



g(«) = 1 

7T 



10 SQ/3 sin(/3G) 



where, 



[10 2sQ + 2 • 10 s " cos(avr) + If / 2 
sin(a7r) 

B = arctan 



(22) 



10 s " + cos(avr) 

and s = log (re /x) with x being the most probable spectral 
parameter. 
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